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1 .  INTRODUCTION  AND  SUMMARY 

The  sensitivity  of  laminates  to  low  velocity  impact  damage  is  a 
practical  design  problem  that  the  ASTM  and  others  have  been  investigating  for 
several  years.  One  of  the  difficulties  in  studying  impact  response  is  predict¬ 
ing  the  complex  interlaminar  stress  gradients  that  it  produces  in  composite 
structures.  This  report  presents  a  systematic  approach  to  interlaminar  stress 
gradient  modeling  with  applications  to  a  graphite-epoxy  laminate  impacted  by  a 
steel  sphere  and  a  laminate  with  an  internal  disbond.  The  approach  is  based 
on  a  family  of  variable  property  finite  elements  for  laminate  models  that 
transition  to  discrete  ply  models  for  regions  with  interlaminar  stress  gradi¬ 
ents.  All  45  elastic  constants  for  a  general  anisotropic  laminate  are  simu¬ 
lated  and  considerable  error  is  shown  to  occur  when  uniform  properties  are 
used  with  certain  laminate  models.  Careful  modeling  of  the  force-deformation 
behavior  is  required  to  predict  accurate  boundary  conditions  in  the  impact 
region  when  bending  is  present.  Numerical  results  are  presented  for  a  thin  8 
ply  laminate  and  a  thick  48  ply  laminate. 

The  report  first  develops  two  new  finite  element  modeling  tools  for 
interlaminar  stress  analysis  and  validates  them  using  control  problems  with 
known  solutions.  The  new  analysis  tools  are  composite  property  models  that 
vary  through  the  thickness  to  correctly  model  bending  and  linear  element 
constraints  for  lamination  theory  force-deformation  modeling.  These  tools  are 
then  used  to  analyze  a  thin  laminate  with  an  internal  disbond  and  a  thick 
laminate  impacted  at  low  velocity  by  a  steel  sphere.  The  disbond  produces  a 
locally  unbalanced  laminate  that  under  biaxial  compression  warps  badly  with 
large  out-of-plane  displacements,  but  with  very  little  internal  force  redis¬ 
tribution.  These  results  were  obtained  from  a  model  using  two  planes  of 
symmetry  that  were  found  to  introduce  ficticious  constraint  forces  between  0 
and  45  degree  plies  at  the  laminate  center.  Results  for  the  laminate  impact 
problem  were  obtained  from  a  complete  360  degree  model  of  the  impact  site.  A 
structural  model  of  the  entire  laminate  was  used  to  obtain  displacement  boun¬ 
dary  conditions  for  a  discrete  model  of  the  impact  site.  These  results  show 
large  subsurface  shears  at  the  perimeter  of  the  contact  area  in  the  second 
ply  below  the  surface.  Only  the  top  four  plies  were  modeled  individually  to 
reduce  costs,  but  more  detailed  models  were  prepared  for  later  solution. 
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2.  INTERLAMINAR  STRESS  FINITE  ELEMENT  MODELING 

Recently  finite  element  modeling  approaches  have  been  suggested  for 
the  determination  of  interlaminar  stresses  in  those  problems  that  cannot  be 
solved  using  lamination  theory  [1,2,3].  All  of  these  involve  some  form  of 
ply-by-ply  modeling  and  an  appeal  to  St.  Venants  principle  to  allow  a  transi¬ 
tion  to  composite  laminate  modeling  away  from  the  region  of  interlaminar 
stress  gradients.  The  issues  affecting  the  accuracy  and  efficiency  of  com¬ 
putational  models  for  this  class  of  problems  are  reviewed  in  this  section  and 
a  sysematic  approach  is  presented  for  use  with  PATCHES-III.  First  a  variable 
property  simulation  of  laminate  force-deformation  behavior  is  introduced  that 
can  represent  all  45  elastic  constants  necessary  to  characterize  general  anis¬ 
otropic  laminate  deformations  [4].  Next  a  family  of  linear  constraint  finite 
elements  are  developed  from  the  PATCHES-III  tricubic  finite  element  to  model 
the  deformations  of  general  laminates.  The  final  and  critical  development 
for  interlaminar  stress  gradient  analysis  is  a  procedure  to  transition  from 
discrete  3D  ply  modeling  to  2D+  composite  laminate  modeling.  The  features  of 
this  new  system  are  then  tested  using  the  edge  effect  demonstration  problem 
and  a  collection  of  balanced  and  unbalanced  laminate  problems. 

2.1  3D  COMPOSITE  PROPERTIES  FOR  LAMINATE  FORCE- 

DEFORMATION  BEHAVIOR 

Pagano  [3]  has  shown  that  the  forces  and  moments  in  a  general  aniso¬ 
tropic  laminate  are  related  to  laminate  strains  and  curvatures  by  45  elastic 
constants  in  the  equations 


*  *  _ 
a.  =  C. .  e.  +  B.  k 
i  lj  j  la  a 


V”2  *  bj 


J'B  j 


*  _ 
e .  +  D 


„  K 
Ba  a 


1  £  1,  j  £  6 


a, 3  =  1 ,2,4 


(2.1) 


which  the  PATCHES-III  contracted  convention  is  used  with  defined  as  the 
in-plane  shear  strain.  These  equations  reduce  to  the  classical  plate  theory 
for  laminates  with  monoclinic  plies  when  =  0.  Once  the  volume  average 
strains  e.-  and  laminate  curvatures  k  have  been  found,  the  individual  ply 
stresses  can  be  determined  using  ply  properties  C^.  and  the  relations  devel¬ 
oped  by  Pagano  [4]. 


HttCSOUO  Him  BUML-NOf  Y1LMD 
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The  elastic  constants  in  Equation  (2.1)  are  functions  of  the 
integrated  moments  of  the  ply  properties 

r  ,  h/2  r  ?  l 

[(A0)ij.»(Al).j.^A2)ij]  =  -fh/2  [1,Z3,Z3J  Cij^z3^dz3  ^2‘2^ 

where  the  exact  expressions  for  che  C. .,  B.  ,  D  0  involve  transformations 

ij  lot  Otp 

arising  from  the  assumption  a3,  c^,  and  are  each  constant  through  the  thick¬ 
ness.  The  approach  developed  in  the  present  study  uses  these  same  moments  to 
define  a  variable  property  composite  material  in  the  thickness  direction  with 

the  same  force-deformation  behavior.  When  a  finite  element  of  this  material 

* 

is  subject  to  the  same  volume  average  strains,  ,  and  curvatures  <a  as  a 
laminate  modeled  using  Equation  (2.1),  the  same  stress  resultants  will  be  pro¬ 
duced.  The  converse  of  this  will  not  be  true  in  general  unless  displacement 
constraints  are  introduced  in  the  thickness  direction  to  prevent  a  pseudo  edge 
effect  from  warping  the  cross-section.  In  this  study,  a  linear  displacement 
constraint  in  the  thickness  direction  is  used  with  the  variable  property  model 
to  define  a  CCL  finite  element  that  predicts  the  force-deformation  behavior  of 
Equation  (2.1). 

Consider  a  parametric  cubic  model  for  the  property  variation  in 
the  normal  direction.  This  allows  up  to  a  sixth  degree  polynomial  in  Z3  when 
the  thickness  coordinate  function,  Z3(£),  is  also  cubic.  However,  the  use  of 
nonuniform  geometry  models  to  aid  in  composite  property  modeling  while  feas¬ 
ible  [Ref.  5]  is  not  very  convenient  as  it  requires  solving  nonlinear  equa¬ 
tions.  The  use  of  a  uniform  (i.e.,  linear)  geometry  model  for  Z3(£), 

Z3(C)  =  h(£  -  1/2)  (2.3) 

and  a  parametric  cubic  for  each  elastic  constant  in  the  thickness  direction, 

c*i (O  =  (S)  C4""1  for  1  s  m  i  4  (2.4) 

.j  ij 

is  adequate  for  most  laminates.  In  some  extreme  cases,  it  may  be  necessary  to 
use  nonuniform  geometry  modeling  to  avoid  unreasonable  values  of  C..j(£;)  at 
local  points  through  the  thickness.  However,  it  is  always  possible  using 
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uniform  models  to  match  the  integrated  moments  in  Equation  (2.2)  which  is 
the  essential  requirement  for  finite  element  modeling.  Substituting  this 
property  model  into  Equation  (2.2)  results  in  four  linear  equations 

(A  )  =  hn+1/(C-l/2)n  C4"mdC(S  ) 
n  ij  0  m  ij 


for  0^ns3 
1  £  mi  4 


These  can  be  expressed  in  matrix  form 

^n  ~  ^nm  3m 

for  any  component  of  C..(^)  where 

*  J 


h/4 

,  h/3  , 

,  h/2  , 

h 

3h2/40 

,  h2/12  , 

,  h2/l 2  , 

0 

7h3/240 

,  h3/30  , 

,  h3/24  , 

h3/l  2 

13h4/l 120 

O 

00 

-C 

* 

,  h4/80  , 

0 

(2.5) 


(2.6) 


(2.7) 


2 

In  the  special  case  of  a  homogeneous  laminate,  A-|  =  A^  =  0  with  =  AqIi  /12, 
the  solution  of  Equation  (2.6)  results  in  simply  C?.(£)  -  (An)../h  as  it  must. 
It  is  possible  to  simulate  Ag.A^^  with  only  a  quadratic,  but  the  use  of  the 
full  cubic  reduces  the  number  of  pathological  cases  for  which  C^(£)  may  not 
be  positive  definite  at  all  points  through  the  thickness.  It  is  important  to 
note  that  this  is  strictly  an  interpolation  problem  that  can  be  avoided  at  the 
expense  of  using  more  interpolation  functions. 

At  this  point,  there  are  two  issues  to  be  resolved  before  the 
approach  can  be  used.  First,  does  the  model  accurately  represent  laminate 
force-deformation  behavior,  and  secondly  is  there  any  significant  difference 
with  simple  constant  property  models  based  on  the  rule  of  mixtures?  The  un¬ 
balanced  laminate  analyzed  by  Pagano  [4]  for  uniaxial  loading  is  used  to 
answer  both  questions  and  to  illustrate  the  practical  limits  of  a  uniform  PC 
property  model . 
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Consider  the  three-ply  laminate  [60/0/-60]  with  ply  properties 

Cn  =  210.35  GPa  Cn  =  C33  =  22.30  GPa 

C12  =  C13  =  7’01  GPa  C23  =  5.75  GPa 

C44  =  C55  =  10.34  GPa  Cg6  =  4.14  GPa 

as  taken  from  Pagano  [2.4]  where  the  msi  units  of  that  paper  are  retained  for 

comparison.  These  ply  properties  result  in  the  force-deformation  properties 

shown  in  Table  1  and  their  inverse  shown  in  Table  2.  Under  uniaxial  load,  for 

example,  the  laminate  twists  and  stretches  with  the  deformations  given  by  the 

corresponding  column  in  Table  2.  Consider  now  the  representation  of  this 

behavior  using  C*.(£)  determined  from  Equation  (2.6).  These  property  functions 
^  3 

in  PC  point  format  are  displayed  in  Table  3.  Note  that  the  extensional  moduli 
were  not  determined  from  Equation  (2.6),  but  are  the  same  constants  as  in 
Table  1.  The  reason  for  this  is  evident  in  Figure  1  which  shows  the  PC  func¬ 
tion  for  C31 (C)-  The  discrete  ply  C^  properties  vary  so  sharply  that  their 
standard  deviation,  14.59,  is  larger  than  their  mean  value,  13.66.  To  repre¬ 
sent  Aq,  A.| ,  and  in  this  case,  a  uniform  PC  function  must  overshoot  zero 
near  the  upper  and  lower  surfaces  which  if  used  would  violate  positive  definite 
C.j j  requirements  at  these  points.  To  avoid  this,  only  AQ  and  A1  were  simu¬ 
lated  by  the  property  model  in  Table  3,  and  as  a  result  the  bending-curvature 
properties  (M-j ,  loading)  are  not  correct.  All  other  force-deformation 
properties  should  be  exact,  and  in  particular  the  uniaxial  case  analyzed  by 
Pagano,  o*  =  1 ,  should  be  modeled  correctly.  This  was  tested  using  one  CCL 
element  in  PATCHES- 1 1 1  with  Table  3  material  properties  and  loaded  by  a  uni¬ 
form  pressure  on  opposite  faces.  The  computed  laminate  strains,  e*  and 
curvatures  k-j  ,  <2>  K4  were  the  same  as  the  Pagano  results.  Table  4.  The 
results  from  a  unit  moment,  M,  =  1.0,  case  are  also  compared  in  Table  4,  and 
these  deformations  are  much  too  low.  The  bending  errors  are  a  direct  conse¬ 
quence  of  the  bending  stiffness  error  caused  by  using  constant  C^  and  C22 
properties  from  the  rule-of -mixtures.  In  this  laminate,  [60/0/-60],  the 
stiffness  error  is  over  100  percent,  and  even  in  a  balanced  structural  lami¬ 
nate,  errors  of  10  percent  to  20  percent  are  common. 

Consider  as  an  example  of  this  behavior  the  48  ply  laminate  analyzed 
later  in  this  report.  The  property  distribution  (c.f.  page  34)  shows  variable 
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Cl 4  and  C24  coupling  through  the  thickness.  This  is  necessary  to  account  for 
the  bending-shearing  coupling  (A2)^  and  ^2)24  in  this  laminate  which  is  quite 
strong.  Note  that  representing  this  behavior  is  important  to  correctly  predict 
the  3D  deformation  response  of  the  laminate  in  areas  of  high  local  bending. 


CnU) 


- c, 


Figure  1.  PC  Property  Modeling  Limits 
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TABLE  1.  THREE-PLY  LAMINATE  FORCE-DEFORMATION  PROPERTIES 

(106  PSI) 


of 

13.656  4. 

,232 

0.925 

-0.694 

°2 

13. 

.656 

0.925 

-1.930 

* 

0„ 

3.234 

-0.018 

4.712  -0.689  -1.925 

0.857 


(SYM)  0.857 


M 

M 

M 


1 

* 

2 

* 

4 


0.494  0.450 

1.500 

0.512 


TABLE  2.  THREE-PLY  LAMINATE  DEFORMATION-FORCE  PROPERTIES 


*  * 

°1  °2 

* 

°3 

★ 

°4 

★ 

0.082  -.018 

-.018 

★ 

e2 

0.167 

-.039 

* 

e3 

0.325 

* 

g4 

0.453 

★ 

g5 

★ 

(SYM) 

* 

Ki 

★ 

k2 

* 

k4 

★ 


°5 


1.167 


qg  m*  Mg  M4 

0.042 

0.602 

0.160 

0.141  0.439 


1.167 

2.829  -.668 
1.559 


4.272 
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TABLE  3.  PARAMETRIC  CUBIC  LAMINATE  PROPERTY  MODEL* 


'11 

'12 

'13 

'14 

'22 

r* 

'23 

r» 

u24 

C33 

C34 

C44 

C55 

C56 

C66 


cV°> 

C^O/3) 

13.658 

13.658 

4.231 

4.231 

.926 

.926 

3.084 

1.827 

13.658 

13.658 

.926 

.926 

8.581 

5.085 

3.234 

3.234 

.078 

.046 

8.284 

3.523 

.550 

1.217 

.385 

.228 

1.550 

.833 

:^<2/3) 

cV> 

13.658 

13.658 

4.231 

4.231 

.926 

.926 

-1.827 

-3.084 

13.658 

13.658 

.926 

.926 

-5.085 

-8.581 

3.234 

3.234 

-.046 

-.078 

3.523 

8.284 

1.217 

.550 

-.228 

-.385 

.833 

1.550 

*Only  extensional  force -deformation  properties  simulated. 


TABLE  4.  FORCE-DEFORMATION  COMPARISONS 


* 

el 

e2 

e* 

e3 

K1 

<2 

k4 

* 

Exact  (o-j  =  1 ) 

.0823 

-.0184 

-.0184 

.0 

.0 

.0 

.0416 

CCL  (a*  =  1) 

.0822 

-.0184 

-.0180 

.0 

.0 

.0 

.0415 

Exact  (M*  =  1) 

.0 

.0 

.0 

.141 

2.829 

-.668 

.0 

CCL  (M*  =  1) 

.0 

.0 

.0 

.065 

.978 

-.193 

.0 

♦Only  extensional  force-deformation  properties  simulated. 
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2.2  CONSTRAINT  FINITE  ELEMENTS 

The  CCL  finite  element  used  to  model  laminate  force-deformation 
behavior  is  one  of  a  family  of  linear  constraint  options  developed  for  the 
PATCHES-III  program.  The  need  for  low  cost  modeling  in  regions  of  uniaxial 
or  biaxial  strain  was  noted  in  an  earlier  report  [6]  and  the  new  elements 
in  Table  5  provide  this  capability.  They  are  based  on  the  linear  constraints 
defined  in  Equation  (2.8) 


P.  =  T.  P 
i  la  a 


where  the  coefficients  T.  are  simply 

la 


1 


T.  = 

i 


0 


i  =  1,2, 3, 4 
a  «  1,4 


2/3  ,  1/3 
1/3  ,  2/3 


1 


(2.8) 


(2.9) 


The  same  coefficients  apply  to  all  three  parametric  coordinates  and,  in  general. 


P .  .  =  T  T  T  P 

ijk  ia  j3  ky  a$y 

If  constraints  are  introduced  in  only  two  coordinates, 

P  =  t  j  £  P 

ijk  ia  jB  k  l  agJt 

and  for  only  one  constraint 

P  =  T  a  £  P 

ijk  ia  js,  km  a£m 


(2.10) 


(2.11) 


(2.12) 


TABLE  5.  PATCHES-III  FINITE  ELEMENT  ADDITIONS 


Displacements* 

Nodes 

Geometry 

Properties 

LLL 

8 

CCC 

CCC 

LLC 

16 

CCC 

CCC 

LCC 

32 

CCC 

CCC 

CCC 

64 

CCC 

ft 

CCC 

*Any  combination  of  L  and  C 

is  avail abl 

e;  L  =  Linear 

C  =  Cubic. 
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Two  key  issues  affecting  the  development  of  the  new  family  of 
elements  are  how  to  efficiently  generate  their  stiffness  matrices  and  how 
to  connect  them  to  each  other.  After  some  early  confusion,  it  was  deter¬ 
mined  that  all  linear  constraints  can  be  applied  before  integration  with  the 
same  result  as  when  they  are  applied  after  integration.  This  greatly  reduces 
the  cost  of  generating  their  stiffness  matrices.  It  is  simple  to  demonstrate 
this  equivalence  in  one  dimension  where  obviously 


Kali 


=  F1<e>  Fj(5)  V* 

■  T1c  FJ<5>  «  fj6 


(2.13) 


but  in  higher  dimensions  interpolatory  quadrature  is  used  in  PATCHES-III  and 
this  caused  some  concern  at  first. 

The  second  issue  was  resolved  by  automating  the  generation  of 
interface  constraints  between  elements  of  different  dimension.  This  allows 
the  user  of  PATCHES-III  to  specify  linear  constraints  on  any  element  or  group 
of  elements  by  simply  placing  a  mnemonic  of  the  type  listed  in  Table  5  on  the 
connectivity  card  for  that  element.  The  program  first  generates  all  explicit 
mesh  point  constraints  and  then  on  a  second  pass  generates  all  interface  or 
implicit  constraints.  This  requires  extensive  checking  for  conflicts,  and  in 
order  to  reduce  their  incidence,  all  three  displacement  components  are  con¬ 
strained  alike.  At  every  constrained  mesh  point,  one  of  the  Equations  (2.10)- 
(2.12)  is  automatically  generated  and  applied  to  all  affected  matrices. 

2.3  INTERLAMINAR  MODELING 

The  use  of  substructuring  to  transition  from  discrete  ply  molding 
to  composite  laminate  modeling  was  recently  Investigated  by  Wang  and  Crossman 
[2].  The  same  technique  has  been  used  for  other  composites  under  the  name 
"minimechanics"  [7]  and,  of  course,  it  has  been  used  for  years  in  aircraft 
structural  analysis.  Substructuring  in  the  case  of  laminates  must  account  for 
two  transition  conditions:  one  along  a  plane  defined  by  a  ply,  and  the  other 
along  a  plane  normal  to  the  laminate.  The  first  case  requires  no  transition 
in  the  shape  of  the  finite  element  mesh  and  appears  to  work  well  for  edge 
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effects  caused  by  uniaxial  loading  [2].  Substructuring  for  the  second  case 
will  require  mesh  changes,  as  well  as  the  transition  to  composite  properties. 

An  approach  using  the  new  CCL  finite  element  with  laminate  properties  is 
evaluated  in  this  section. 

Consider  the  interlaminar  stress  problem  analyzed  in  the  PATCHES-III 
User's  Manual.  First,  the  accuracy  of  the  new  constraint  finite  elements  will 
be  demonstrated  for  the  same  four  element  model  using  the  original  properties. 
The  two  interior  elements  were  constrained  to  be  all  linear  (LLL)  and  the  two 
edge  elements  were  constrained  linear  in  the  direction  of  the  load  (LCC).  No 
significant  change  in  the  earlier  results  should  occur  because  the  displacement 
solution  has  the  same  form  in  these  regions.  Normal  stress  comparisons  in 
Table  6  show  this  to  be  true  even  though  the  model  has  been  reduced  to  only  92 
degrees-of-freedom.  In  fact,  the  assumption  of  linear  displacements  in  the 
thickness  direction  (LCL)  also  was  used  with  very  little  mid-surface  normal 
stress  error  as  comparisons  in  Figure  2  show.  Interlaminar  stress  gradients 
between  the  0°  and  90°  plies,  however,  are  drastically  changed  in  the  LCL  model 
as  the  shear  stress  comparisons  in  Table  7  show.  Note  that  the  interlaminar 
shear  stress  appears  to  contain  a  singularity  at  the  free-edge  between  the  0° 
and  90°  plies,  023  (o,b,h).  The  uniform  mesh  model  gives  the  mean  value  of 
©23,  but  one  element  cannot  model  this  response  and  return  to  zero  at  the  free- 
edge.  A  nonuniform  mesh  element  does  a  better  job  of  matching  this  condition 
and  produces  about  as  much  accuracy  as  can  be  expected  from  a  single  LCC 
element. 


TABLE  6. 

INTERLAMINAR 

NORMAL  STRESS 

COMPARISONS* 

★  * 

C/2h 

ccc/Ccc 

LLL/LCC 

LLL/LCL 

(428  D.O.F.)  (92  D.O.F.) 

(44  D.O.F.) 

0 

+2.95 

+2.89 

+2.76 

1/3 

-  .26 

-  .27 

-  .31 

2/3 

-  .43 

-  .44 

-  .25 

1 

-  .16 

-  .05 

-  .03 

2 

-  .02 

-  .03 

-  .03 

3 

+  .01 

-  .01 

-  .00 

4 

-  .01 

+  .02 

-  .02 

Comparisons  at  midsurface  between  90°  plies. 
**Z2  =  b  -  distance  from  free-edge. 
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TABLE  7.  INTERLAMINAR  SHEAR  STRESS  COMPARISONS* 


£/2h 

CCC/CCC 

Uniform 

LLL/LCC 

Uniform 

LLL/LCL 

Uniform 

C/2h 

LLL/LCC 

Nonuniform** 

0 

-2.72 

-2.72 

-.33 

0 

+  .58 

1/3 

-  .85 

-  .85 

-.47 

1/9 

-1.94 

2/3 

-  .60 

-  .60 

-.42 

4/9 

-  .67 

1 

+  .03 

+  .03 

+  .03 

1 

-  .31 

2 

-  .05 

-  .05 

-.04 

2 

-  .03 

3 

-  .02 

-  .02 

-.02 

3 

-  .02 

4 

.00 

.00 

.00 

4 

.00 

♦Comparisons  between  90°  and  0°  plies. 
**Z£  =  6  +  4£  -  (edge  elements  only). 


Consider  now  the  transition  from  discrete  ply  properties  to  a 
composite  property  model  C*j  away  from  the  free-edge.  A  six  element  model 
shown  in  Figure  3  was  used  to  determine  model  accuracy  when  transitions  of 
this  type  are  made  using  the  present  constraint  finite  elements.  Note  the 
two  edge  elements  are  exactly  as  before,  and  the  transition  mesh  uses  three 
wedge  shaped  elements.  Rule-of-mixture  composite  properties  were  used  for 
elements  number  one  and  three,  and  all  other  elements  use  local  ply  properties 
as  in  the  original  model.  The  interlaminar  stresses  at  the  free  edge  changed 
very  little.  Figure  2.  At  the  interface,  of  course,  small  local  shear  stresses 
occur  because  of  the  discrete  change  in  properties.  Their  maximum  value  is 
less  than  one-half  of  one  percent  the  applied  axial  stress.  Essentially,  the 
same  results  were  obtained  using  C*j(£)  variable  composite  properties.  Inter¬ 
estingly,  the  two  ply  laminate  properties  were  easier  to  model  with  a  cubic 
than  the  three  ply  laminate  because  the  [0/90]  properties  are  an  odd  function. 
These  results  indicate  the  transition  from  discrete  ply  modeling  to  composite 
laminate  modeling  can  be  made  using  the  new  constraint  elements. 


1 

LCC 

LAMINATE  ,  LLC 

2 

LCC 

PLY  ,  90° 

3 

LCC 

LAMINATE  ,  LLC 

4 

LCC 

PLY  ,  0° 

5 

LCC 

PLY  ,  90° 

6 

LCC 

PLY  ,  0° 

Figure  3.  PATCHES- III  Transition  Model  for  Interlaminar 
Stress  Analysis 
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3.  DISBONDED  LAMINATE  ANALYSIS 

The  effects  of  an  unsymmetric  disbond  on  laminate  deformations  and 
associated  internal  load  redistributions  around  the  disbond  are  investigated 
for  a  thin  laminate.  An  important  aspect  of  the  unsymmetrically  disbonded 
laminate  problem  not  considered  in  the  present  investigation  is  the  effect  of 
laminate  thickness.  There  is  considerable  experimental  evidence,  Williams  [8, 
9],  that  the  strength  of  thick  laminates  (-50  plies  or  more)  is  sensitive  to 
internal  disbonds  in  both  unidirectional  tape  and  bidirectional  fabric  mate¬ 
rials  of  graphite-epoxy.  These  remarks  refer  to  compressive  not  tensile  loads 
and  to  disbonds  between  plies  that  occur  without  fiber  damage.  Results  from 
the  present  investigation  suggest  thin  laminates  (~10  plies  or  less)  will  ex¬ 
perience  a  loss  in  stiffness  caused  by  out-of-plane  bending  before  any  stress 
critical  failure.  The  laminate  in  this  study  is  an  eight  ply  graphite-epoxy 
flat  panel  with  a  small  circular  disbond  (T/R  =  0.1)  located  three  plies  from 
the  surface  and  at  the  center  of  the  panel.  A  ply-by-ply  modeling  of  this 
problem  is  used  to  study  interlaminar  force-deformation  behavior  in  the  vicin¬ 
ity  of  the  disbond  with  CCL  finite  elements.  No  attempt  was  made  to  include 
in  the  model  any  strain  singularities  that  might  exist,  although  this  could 
have  been  done  had  fracture  mechanics  been  the  focus  of  the  study.  The  objec¬ 
tive  was  to  determine  the  deformations  and  internal  load  redistribution  caused 
by  an  unsymmetric  disbond  in  a  compressively  loaded  thin  laminate. 

The  disbonded  laminate  problem  was  analyzed  early  in  the  study  using 
symmetry  boundary  conditions  on  two  planes.  Prototype  models  indicated  a 
definite  skewed  displacement  response  would  occur  and  that  the  use  of  symmetry 
planes  would  inhibit  some  response  modes.  Unfortunately,  PATCHES-III  was 
limited  to  50  elements  at  that  time  and  no  practical  means  of  avoiding  symmetry 
modeling  was  available.  The  resulting  symmetry  model  solution  is  essentially 
a  3D  laminate  solution  with  interlaminar  stress  gradients  distorted  at  the 
origin  by  the  symmetry  boundary  conditions.  However,  the  force-deformation 
behavior  of  the  laminate  is  represented  reasonably  well  and  several  interest¬ 
ing  conclusions  can  be  drawn  from  the  symmetry  solution. 

3.1  PROTOTYPE  MODELS 

A  schematic  of  the  Idealized  disbonded  laminate  analyzed  In  this 
study  is  shown  in  Figure  4.  The  disbonded  region  is  at  the  center  of  a  square 
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laminate  between  plies  five  (-45°)  and  six  (0°).  The  region  is  circular  with 
a  bubble  shape  whose  profile  is  given  by 

Z3  =  tp(1+cos  nr/R)/2  (3.1) 

where  tp  is  a  ply  thickness  and  the  maximum  radius  is  R  =  1.27  cm  (0.5  in.). 
Note  that  this  shape  is  tangent  to  the  adjacent  bonded  region  and  corresponds 
to  an  elastically  warped  surface.  A  catenary  shape  would  correspond  to 
inelastic  fiber  damage  at  the  boundary  of  the  disbonded  region.  The  laminate 
is  15.24  cm  (6  in.)  square  and  uniformly  compressed  0.01524  cm  (.006  in.)  in 
both  the  e^  and  e^  directions.  At  the  loaded  edges,  the  normal  displacement 
component,  u^,  is  restrained  to  zero.  The  disbonded  plies  [-45/45/0]  are 
unbalanced  and  after  their  separation,  the  remaining  plies  [O/45/-45/O2]  are 
also  unbalanced.  As  in  the  unbalanced  three-ply  laminate  analyzed  by  Pagano, 
there  will  be  stretching-twisting  coupling  that  must  be  considered  in  select¬ 
ing  a  symmetry  model.  It  is  also  important  to  examine  the  inplane  force 
changes  caused  by  a  disbond  to  assist  in  preparation  of  the  ply-by-ply  model 
of  the  laminate. 

A  small  two  element  model  was  used  to  examine  interply  force  changes 
when  a  -45°  ply  and  a  0°  ply  disbond  locally.  The  prototype  model.  Figure  5, 
consists  of  a  -45°  ply  and  a  0°  ply  in  one  quadrant  of  the  laminate  (7.62  cm) 
loaded  by  a  uniform  strain  =  0.001.  To  simulate  the  disbond,  the  two  ele¬ 
ments  were  disconnected  at  node  8  corresponding  to  the  center  of  the  laminate. 
The  two  elements  were  type  LLL  and  the  nodal  forces  at  the  interply  corners 
were  used  to  measure  the  effect  of  the  disbond.  These  are  shown  in  Figure  6 
and  indicate  only  a  small  constraint  force  is  relaxed  by  a  disbond  at  node  8. 
Note,  however,  that  the  effect  is  to  increase  the  transverse  load  in  the  0°  ply 
which  carries  most  of  the  load  associated  with  the  applied  strain.  It  is 
obvious  that  a  larger  change  will  occur  if  either  nodes  5  or  7  are  cut  instead 
of  node  8  which  is  analogous  to  a  disbond  at  node  8  in  a  +45/0°  layup.  Results 
from  this  case.  Figure  7,  do  show  a  larger  change,  but  it  reduces  the  trans¬ 
verse  load  in  the  0°  ply.  The  directional  dependence  in  stiffness  produces  a 
directional  dependence  in  the  effect  of  a  disbond  under  any  given  load.  In 
this  case  of  a  circular  disbond,  Figure  4,  under  biaxial  load,  it  appears 
that  interply  constraint  forces  will  vary  around  the  perimeter  with  relative 
maximums  or  minimums  every  45°. 
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El  =  137.90  GPa 
Et  =  9.31  GPa 
G  =  9.31  GPa 
V  =  0.3 


8,13 

Disbonded  Case 


u  •  e-j  =  0  on  Faces  4-8-5-1 ,  13-12-9-5 
u  •  e-j  =  0.03  on  Faces  3-7-6-2,  7-11-10-6 


1 


Figure  5.  Two-Ply  Disbond  Prototype  Model  for  Forces 
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Next,  a  three-ply  [0/45/-45]  laminate  model  was  used  to  examine  out- 
of-plane  warping  under  inplane  loading  of  the  unbalanced  disbond  material. 

Each  ply  was  modeled  with  CCL  elements,  and  the  dimension  of  the  square  model 
was  made  equal  the  radius  of  the  disbond.  Under  uniaxial  load  there  is  elastic 
coupling  between  stretching  and  twisting  that  caused  large  out-of-plane  dis¬ 
placements  because  of  the  low  bending  stiffness  of  a  three-ply  laminate  (the 
plies  are  only  0.01778  cm  thick).  The  stresses,  strains  and  deformations  of 
this  [0/45/-45]  laminate  are  all  either  asymmetric  or  symmetric  about  axes  at 
±45°  to  axes  e-| ,  e£.  Figure  8.  These  results  and  the  [0/45]  constraint  force 
results  indicate  diagonal  symmetry  rather  than  Cartesian  symmetry  should  be 
used  in  modeling  force-deformation  behavior  of  the  disbonded  laminate. 

Another  observation  from  the  unbalanced  prototype  models  was  the 
extremely  slow  convergence  of  the  conjugate-gradient  solution  procedure.  The 
membrane  response  accounts  for  99.7  percent  of  the  potential  energy  and  this 
result  is  converged  in  less  than  N  cycles.  However,  the  low  energy  out-of¬ 
plane  bending  response  required  almost  4N  cycles  to  converge  with  little  change 
to  principal  stresses  or  strains.  Note  that  the  interply  normal  forces  are  two 
orders  of  magnitude  smaller  than  the  inplane  forces.  The  models  were  changed 
to  CCC  elements  to  check  for  purely  numerical  ill -conditioning  and  the  same 
slow  convergence  was  observed.  No  solution  to  the  pathological  computational 
behavior  caused  by  weak  out-of-plane  coupling  was  found. 

3.2  LAMINATE  FINITE  ELEMENT  MODELS 

The  dimensions  of  the  laminate  were  specified  to  insure  no  coupling 
between  the  displacement  boundary  conditions  and  deformations  in  the  vicinity 
of  the  disbond.  The  laminate  is  square  of  width  15.24  cm  (6  inches)  with  a 
disbond  radius  of  1.27  cm  (.5-inch)  between  plies  five  and  six.  Figure  4.  The 
plies  are  graphite-epoxy  of  thickness  0.01778  cm  (0.007-inch)  with  their 
nominal  elastic  properties  as  given  in  Table  8.  The  stacking  sequence  for  the 
laminate  is  [0,45,-45,02,-45,45,0]j  with  the  disbond  between  a  0°  ply  and  a 
-45°  ply.  The  edges  are  loaded  by  imposing  a  uniform  inplane  displacement  of 
0.00762  cm  (0.003-inch)  which  produces  biaxial  compression. 

A  control  model  of  the  laminate  was  constructed  using  the  same  mesh 
in  the  plane  as  in  Figure  4,  but  with  only  four  plies  [0,45,-45,0]^ 

through  the  half  thickness.  The  purpose  of  this  model  was  to  provide  before- 
and-after  data  for  comparison  with,  the  disbonded  results,  in  particular  the 
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Figure  8.  Unbalanced  Laminate  Bending  Displacements 
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TABLE  8.  GRAPHITE-EPOXY  NOMINAL  ELASTIC  PROPERTIES 


En  =  137.90 

GPa  (20  msi) 

v12  =  0.3 

G12 

=4.83  GPa  (.7  msi ) 

E22  =  9.31 

GPa  (1.35  msi ) 

V1 3  =  0,3 

G1 3 

=  4.83  GPa  (.7  msi) 

E33  =  9'31 

GPa  (1 .35  msi) 

v23  =  0,3 

G23 

=  4.83  GPa  (.7  msi) 

internal  forces.  At  the  same  time  the  control  model  shows  the  magnitude  of 
the  interlaminar  constraint  forces  caused  by  the  use  of  diagonal  symmetry, 

Table  9.  The  forces  F 2  should  be  exactly  zero  if  the  interlaminar  deformations 
are  zero  at  the  center  of  the  laminate.  Instead,  a  large  self-equilibrating 
pair  of  forces  occur  between  the  0°  plies  and  the  ±45°  plies.  These  forces 
distort  the  interlaminar  stresses  at  the  origin.  It  should  be  noted  that  F 2 
sums  to  zero  between  the  +45  and  -45  and  the  use  of  diagonal  symmetry  for  these 
plies  is  consistent  with  material  symmetry  axes. 

TABLE  9.  BONDED  LAMINATE  CONSTRAINT  FORCES  THROUGH 
THE  THICKNESS  AT  r  =  0 


Z,  =  4tn 

3  p 

Zn  =  3tn 

3  p 

U  =  2tn 

3  p 

Z,  =  tn 

3  p 

Z-.  = 

3  p 

NODE 

49 

50 

51 

52 

53 

F1 

54.98  N 

100.44  N 

102.00  N 

100.22  N 

50.18  N 

f2 

0. 

43.41  N 

0. 

-43.41  N 

0. 

The  control  model  also  demonstrates  that  even  the  small  biaxial 
strain  (e  s  0.1%)  imposed  at  the  edges  of  the  panel  is  sufficient  to  buckle 
the  laminate  elastically.  An  engineering  check  was  made  using  isotropic  sim¬ 
ply  supported  plate  formulas  with  the  stiff  direction  properties.  This  cal¬ 
culation  indicates  the  laminate  will  buckle  at  P^  =  P2  ~  2000  N,  while  the 
applied  load  is  P  >  10,000  N.  The  exact  solution  for  buckling  of  a  laminated 

anisotropic  plate  [10]  would  be  lower  because  the  small  number  of  plies  accen¬ 
tuates  the  bending-stretching  coupling.  It  seems  likely  that  the  disbond  will 

reduce  the  buckling  load  even  further  or  possibly  cause  a  large  deflection 

stability  problem  in  a  laminate  this  thin,  T  =  0.142  cm  (0.056-inch). 
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3.3  DISBONDED  LAMINATE  RESULTS 

The  effect  of  the  disbond  on  the  interlaminar  force  distribution  is 
quite  small  and  the  effect  on  bending  deformation  is  quite  large.  The  center- 
point  initial  displacement  caused  by  the  disbond  is  0.5  tp  before  any  load  is 
applied.  The  centerpoint  elastic  displacement  after  loading  to  a  biaxial 
strain  of  only  0.1  percent  is  an  additional  0.4  tp.  Figure  9.  In  other  words, 
a  unit  imposed  displacement  (biaxial)  in  the  plane  of  the  laminate  produces 
almost  a  unit  normal  displacement  at  the  center  of  the  laminate.  This  result 
suggests  the  disbond  will  initially  lead  to  a  stiffness  critical  failure 
rather  than  a  stress  failure.  The  reason  for  this  behavior  seems  to  be  the 
extremely  low  bending  stiffness  which  produces  large  out-of-plane  deformation 
when  the  disbond  couples  the  membrane  and  bending  displacement  response.  It 
is  important  to  note  that  the  potential  energy  change  between  the  bonded  and 
disbonded  cases  is  less  that  one  milli  Joule,  Table  10,  which  corresponds  to 
an  energy  release  rate  of  less  than  1  N/m.  In  contrast,  the  critical  energy 
release  rate  for  an  edge  delamination.  Reference  11,  is  will  over  100  N/m  for 
an  11  ply  graphite-epoxy  laminate. 

TABLE  10.  DISBONDED  LAMINATE  ENERGY  CHANGE 


Potential 

Energy 

Bonded 

Disbonded 

3.958455  Joules 

3.958069  Joules 

(35.03528 

(35.03186 

in. -lbs) 
in.-ltjs) 

The  small  redistribution  of  internal  forces  is  evident  in  Table  11 
which  shows  the  center  point  constraint  forces  through  the  thickness.  Note 
that  the  forces  at  nodes  53  and  54  are  summed  when  the  disbond  is  closed. 

Note  also  that  the  force  at  node  53  in  the  bonded  model  must  be  added  to  node 
52  for  comparison  with  the  disbond  results.  This  is  because  only  one  element 
was  used  for  the  two  center  plies.  Figure  4,  in  the  disbond  model.  A  check 
of  the  force  changes  through  the  thickness  was  also  made  at  r  =  R,  Table  12. 

To  a  large  extent  the  symmetry  boundary  conditions  inhibit  in-plane 
redistribution  and  to  some  extent  interlaminar  normal  force  changes.  Con¬ 
sidering  the  large  normal  displacements.  Figure  9,  one  might  expect  an 
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Figure  9.  Normal  Deformations  in  the  Plies  Adjacent  to  the  Disbond 
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TABLE  11.  DISBONDED  LAMINATE  CONSTRAINT  FORCE  CHANGES 
THROUGH  THE  THICKNESS  AT  r  =  0 


Node 

Bonded 

Disbonded 

Z3 

F1 

F2 

Z3 

F1 

F2 

49 

4(p 

54.98  N 

0.  N 

4tP 

44.54  N 

1.01  N 

50 

3tP 

100.44  N 

43.41  N 

*p 

103.68  N 

45.75  N 

51 

102.00  N 

0.  N 

105.69  N 

-.22  N 

52 

*p 

100.22  N 

-43.41  N 

*P 

160.78  N 

-46.58  N 

53* 

0 

50.18  N 

0.  N 

-V 

112.41  N 

-.90  N 

54 

-1.5  tp 

43.62  N 

-41.33  N 

55 

(Symmetric) 

-2- 5 

90.67  N 

3.22  N 

56 

-3.5  tp 

89.49  N 

38.57  N 

57 

-4.5  tn 

42.46 

1.43  N 

P 

*Node  53  is  at  Zj  =  0  in  the  bonded  laminate  case. 


interlaminar  normal  force  to  occur  at  r  «  R.  As  the  data  in  Table  12  show, 
only  a  small  increase  in  the  small  interlaminar  normal  force  occurred  at  node 
41.  A  possible  explanation  is  that  the  symmetry  boundary  conditions  at  the 
center  prevents  nodes  53  and  54  from  moving  relative  to  each  other  in-plane 
as  they  probably  would  in  a  complete  solution.  If  such  deformations  occur, 
they  would  create  transverse  shear  gradients  in-plane  and  equilibrium  would 
require  a  normal  stress  gradient  in  the  thickness  direction, 

°31,1  +  °32,2  +  a33,3  =  0 

The  magnitude  of  the  spurious  force  at  node  54  suggests  this  deformation 
could  be  significant  and  would  certainly  increase  the  energy  release  rate. 
However,  it  would  have  to  increase  by  a  factor  of  100  before  the  delamination 
would  propagate  if  the  data  in  Reference  11  apply.  This  uncertainty  can  only 
be  resolved  by  analyzing  a  complete  model  without  the  assumption  of  symmetry 
conditions. 
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TABLE  12.  DISBONDED  LAMINATE  CONSTRAINT  FORCE  CHANGES 
THROUGH  THE  THICKNESS  AT  r  =  R 


Node 

Bonded 

Disbonded 

Z3 

F1  *  F2 

F  * 
r3 

Z3 

F1  =  F2 

F3 

33 

4S 

12.18  N 

0.  N 

4tP 

14.01  N 

0.  N 

35 

74.54  N 

±0.56  N 

3tP 

75.17  N 

±0.07  N 

37 

68.29  N 

±1.12  N 

*S 

69.10  N 

±0.14  N 

39 

S 

33.09  N 

±0.84  N 

1p 

33.13  N 

±0.93  N 

41 

-*p 

-{p 

32.16  N 

±1.49  N 

43 

-2tP 

(Symmetric) 

-2tP 

67.28  N 

±0.84  N 

45 

-3tP 

-3tp 

73.04  N 

±0.34  N 

47 

-4t 

12.19  N 

0.  N 

P 

P 

*Fg  forces  are  self-equilibrating  interlaminar  forces. 


The  interlaminar  stress  results  in  keeping  with  the  force  results 
show  only  small  changes  in  their  distribution.  On  a  percentage  basis,  the 
largest  change  occurs  in  the  bottom  ply  where  the  small  oLT  shear  stress  is 
doubled.  A  plot  of  this  stress  component  around  the  disbond  at  r  =  R,  Figure 
10,  shows  relative  maximums  and  minimums  at  ±45°  as  expected.  The  use  of 
symmetry  boundary  conditions  also  caused  local  distortions  at  the  origin  in 
the  stress  results  that  make  their  magnitudes  questionable.  However,  these 
conditions  affect  the  bonded  and  disbonded  cases  equally. 


? 

[ 
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Figure  10.  Botto*  Ply  Shear  Stress  Changes 
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4.  LAMINATE  IMPACT  ANALYSES 

At  low  velocities  impact  damage  can  be  modeled  as  a  quasi-static 
response  problem  using  Hertzian  contact  pressures  in  most  cases.  This  approach 
is  used  in  the  present  study,  and  Greszczuk,  Reference  12,  gives  an  excellent 
account  of  the  method  and  its  rationale.  The  only  analytic  results  available 
are  for  the  elastic  half  space  problem  which  was  recently  solved  for  a  trans¬ 
versely  isotropic  material.  Reference  13.  One  effect  of  transverse  isotropy 
shown  graphically  in  that  paper  is  a  shift  of  the  maximum  Von  Mises  stress 
from  the  impact  centerline  to  a  radius  approaching  the  contact  radius.  The 
significance  of  this  stress  parameter  for  a  composite  is  dubious,  but  it  does 
indicate  a  trend  of  increasing  transverse  shear  below  the  contact  radius, 
which  also  occurs  in  the  present  results. 

The  disbonded  laminate  analysis  problems  caused  by  symmetry  displace¬ 
ment  boundary  conditions  led  to  the  use  of  a  full  360  degree  finite  element 
model  for  the  impact  problem.  Another  factor  in  this  decision  was  the  earlier 
PATCHES-III  sandwich  panel  impact  analysis  [6]  which  also  showed  spurious  con¬ 
straint  forces  between  the  0°  and  45°  plies  on  the  centerline  of  a  double  sym¬ 
metry  model.  In  order  to  model  a  full  360  degrees,  ply-by-ply,  the  program 
was  modified  to  allow  several  hundred  solid  elements.  Unfortunately,  the  com¬ 
puter  budget  available  was  sufficient  to  analyze  only  a  32  element  model  after 
all  computational  problems  were  resolved.  These  results,  however,  provide  the 
first  in-depth  picture  of  the  asymmetric  interlaminar  stress  gradients  that 
occur  in  a  thick  laminate  impacted  off-center.  More  detailed  models  were  pre¬ 
pared  for  future  analysis  and  these  are  in  an  appendix  to  this  report. 

4.1  THICK  LAMINATE  COMPOSITE  RESPONSE 

A  48  ply  graphite-epoxy  laminate  supported  between  two  rigid  spars 
is  shown  schematically  in  Figure  11  being  impacted  by  a  steel  sphere.  The 
impact  point  is  close  to  a  spar  and  in  this  problem  a  contact  force  of  4448  N 
was  specified.  Considering  the  large  number  of  plies  in  this  problem,  not  all 
plies  can  be  modeled  for  any  given  analysis.  This  practical  consideration 
becomes  even  more  critical  when  symmetry  conditions  cannot  be  used.  First  a 
variable  property  30  laminate  model  was  used  to  determine  the  bending  deforma¬ 
tions  in  the  vicinity  of  the  impact  site.  This  analysis  used  a  large  but 
finite  width  strip  to  avoid  boundary  condition  effects,  Figure  12.  The  30 
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laminate  analysis  produced  large  transverse  shear  stresses  at  the  impact  site 
approaching  the  bending  stresses  in  magnitude.  Laminate  force-deformation 
properties  were  modeled  using  the  C?.(£)  formulation  derived  earlier.  Table  13. 

*  J 

As  these  data  indicate,  only  the  off-diagonal  coupling  terms  vary  appreciably 
through  the  thickness. 


LAMINATE:  [±45/02/±45/02/±45/0/90]2s 
GRAPHITE-EPOXY 


Figure  11.  Thick  Laminate  Impact  Schematic 
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Laminate  Structural  Response  PATCHES-III  Model 
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TABLE  13.  PARAMETRIC  CUBIC  THICK  LAMINATE  PROPERTY  MODEL* 


*  clj<» 

C* j(l/3)  =  ctj(2/3) 

Cll 

72.622 

71.650 

Cl  2 

21.712 

14.906 

Cl  3 

2.999 

2.910 

C14 

2.427 

-0.807 

C22 

24.173 

36.887 

C23 

2.379 

2.468 

C24 

2.427 

-0.807 

C33 

9.432 

9.432 

C44 

22.091 

16.279 

C55 

3.771 

3.771 

C66 

3.744 

3.744 

♦Stiffness  in  GPa  units  (1  Psi  =  6894.757  Pascals) 


4.2  IMPACT  SITE  MODEL 

The  composite  laminate  displacement  response  for  element  number  one 
Figure  12,  provided  boundary  conditions  for  the  impact  site  model  shown  in 
Figure  13.  The  contact  radius  was  determined  using  both  the  Greszczuk  [12] 
equations  and  the  Dahan  and  Zarka  [13]  equations.  The  latter  have  typograph¬ 
ical  errors  that  were  reconciled  using  the  limiting  case  of  an  isotropic  half 
space  to  yield 


r 


o 


(4.1) 


where  6^  and  62  are  lengthy  algebraic  functions  of  the  half-space  flexibili¬ 
ties  St..  Substituting  the  laminate  St.  in  these  expressions  gave 


61  =  -8.06  x  10“12  m2/N 

62  =  -68.09  x  10"12  m2/N 


The  remaining  terms  in  Equation  (4.1)  are  the  force  F  *  4448  N,  the  sphere 
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Figor®  13.  Impact  Sit®  PATCHES-lIl  M®dd 
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radius  R  =  6.35  mm  and  the  elastic  constants  E  =  206.84  Gpa  and  v  =  0.3  for  the 
sphere.  These  data  produced  a  contact  radius  of  r  =  1.134  mm.  The  Greszczuk 
equation  for  contact  radius 

ro  =  [t  FR  lkl  +  k2>]  V3  <4-2> 

also  uses  the  S*j  to  compute  related  coefficients  and  k2  which  are 

k1  =  1.4004  x  10"12  m2/N 
k2  =  19.8173  x  10-12  m2/N 

producing  a  contact  radius  of  1.122  mm.  The  smaller  radius  was  used  to  model 
the  impact  site  so  that  T/R  =  5.65. 

Several  disastrous  attempts  at  analyzing  the  impact  site  model, 

Figure  13,  were  made  using  CCC  elements  for  the  upper  plies  and  CCL  elements 
for  the  subsurface  laminate  model.  These  attempts  never  quite  converged  and 
the  deformations  were  highly  distorted.  After  considerable  effort  it  was 
determined  that  the  model  had  pathological  stiffness  discontinuities  in  the 
thickness  direction  that  were  easily  removed  once  found.  At  the  impact  site 
the  elements  through  the  thickness  must  be  either  all  CCC  or  all  CCL  to  avoid 
this  condition.  When  the  two  element  types  are  mixed  under  impact  loading 
conditions,  the  CCC  elements  deform  as  if  they  were  pressed  against  a  rigid 
boundary.  After  resolving  this  modeling  dilemma,  excellent  results  were 
obtained  using  all  CCC  elements. 

4.3  INTERLAMINAR  IMPACT  STRESS  RESULTS 

It  is  important  to  keep  in  mind  that  the  present  results  are  for  a 
smal 1  sphere  impacting  a  thick  laminate,  T/R  =  5.65,  with  sufficient  energy  to 
generate  a  contact  force  of  4448  N  (1000  lbs).  The  resulting  stress-strain 
gradients  are  quite  high  in  the  upper  ply  group.  Figure  14,  especially  in  mate¬ 
rial  coordinates.  Note  in  particular  that  the  second  ply,  a  -  45,  has  larger 
fiber  stress  gradients  in  the  thickness  direction  than  the  top  ply,  although 
the  magnitudes  are  slightly  lower.  The  centerline  Cartesian  shear  strains. 
Figure  15,  show  rather  dramatically  that  the  use  of  two  symmetry  planes  would 
have  produced  large  shear  errors  along  the  centerline.  Interesting,  the 
chordwise  shear  strain,  ££3,  is  very  small  indicating  one  plane  of  symmetry 
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Figure  15.  Cartesian  Centerline  Shear  Strains 
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in  the  response.  If  the  loading  were  not  in  the  transverse  direction,  or  if 
the  laminate  were  unbalanced,  even  this  symmetry  condition  would  not  exist. 

Consider  next  the  stress  distribution  around  the  perimeter  of  the 
impact  site.  The  fiber  stresses  at  the  top  surface  as  a  function  of  circum¬ 
ferential  position,  Oii (0) »  are  shown  in  Figure  16.  The  tension  stresses  are 
confined  to  a  thin  layer  at  the  free  surface  and  are  probably  not  a  good 
indicator  of  actual  fiber  stresses.  When  the  fiber  stresses  are  averaged 
through  a  single  ply  thickness.  Figure  16,  the  values  are  all  compressive. 

In  either  case,  the  distribution  is  periodic  in  9  with  a  period  equal  ir.  The 
distribution  of  transverse  shear  stress  in  material  coordinates  is  also  peri¬ 
odic,  Figure  17,  but  with  period  equal  2tt.  The  two  shear  stress  components 
in  this  figure,  o^3  and  ©23,  were  averaged  over  one  ply  thickness,  the  second 
ply.  This  ply  has  the  maximum  transverse  shear  in  the  laminate  and  it  occurs 
at  the  impact  perimeter  not  the  centerline.  Similar  behavior  was  observed  by 
Dahan  and  Zarka  in  zinc  and  cadmium  during  elastic  contact  with  a  steel  sphere. 
These  materials  also  had  their  low  stiffness  in  the  direction  normal  to  the 
impact  surface,  but  they  are  not  as  anisotropic  as  the  laminate  material.  It 
is  also  interesting  to  note  that  the  two  shear  stress  components  0^3  and  ^3 
are  phase  shifted  by  90°  in  the  circumferential  direction.  The  in-plane  ply 
shear  stress  in  Figure  18  shows  a  quasi-periodic  distribution  of  period 
equal  it  that  is  similar  to  the  fiber  stress  distribution.  However,  the  struc¬ 
tural  response  of  the  laminate  appears  to  be  superimposed  on  the  elastic  half¬ 
space  response.  This  would  explain  the  Peak  nearest  the  center  of  the 
laminate  being  slightly  smaller.  There  is  also  a  subharmonic  component  that 
may  be  caused  by  the  composite  solution  boundary  conditions  applied  to  the 
impact  site  model.  To  explore  this  second  order  effect  would  require  addi¬ 
tional  analyses.  A  summary  of  the  ply  material  stress  maximums  is  provided 
in  Table  14. 

It  is  interesting  to  observe  the  similarities  and  differences  in  a 
3D  laminate  solution  and  a  3D  elasticity  solution  in  the  impact  area.  The 
deformations  at  the  top  surface  are  compared  in  Figure  19  and  show  very  simi¬ 
lar  shapes.  The  maximum  3D  laminate  displacement  is  about  85  percent  of  the 
elasticity  result;  however,  the  laminate  strains  are  an  order  of  magnitude  too 
low  at  the  top  surface.  The  elasticity  solution  using  laminate  properties  was 
much  closer  but  still  low.  The  normal  strain  in  this  case  is  quite  close,  but 
the  in-plane  strains  are  one-half  to  one-third  the  ply  strain  results. 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 

A  systematic  approach  to  modeling  complex  interlaminar  stress  gra¬ 
dients  has  been  developed  and  applied  to  two  graphite-epoxy  laminates.  The 
approach  uses  variable  property  3D  finite  elements  that  correctly  model  force- 
deformation  behavior  and  constant  property  3D  finite  elements  that  accurately 
model  ply  shear  strains.  Analyses  for  this  class  of  problems  are  difficult, 
but  considerable  progress  has  been  made  toward  solving  the  computational  and 
modeling  issues  encountered. 

•  To  model  composite  laminate  force-deformation  behavior 
with  3D  finite  elements  requires,  in  general,  proper¬ 
ties  that  vary  through  the  laminate  thickness. 

•  Interlaminar  deformations  are  not  symmetric  in  a 
balanced  symmetric  laminate  under  symmetric  in-plane 
loads. 

•  Interlaminar  deformations  have  one  plane  of  symmetry  in 
a  balanced  symmetric  laminate  under  normal  Hertzian 
contact  pressure. 

Analysis  of  the  disbonded  laminate  force  deformation  behavior  indi¬ 
cates  a  large  out-of-plane  displacement  response  occurs  without  singificant 
internal  force  redistribution.  Comparison  with  a  bonded  laminate  control 
model  also  indicates  very  little  energy  is  released  by  the  disbond.  The  new 
CCL  finite  element  worked  well  in  this  application  and  should  be  an  effective 
new  tool  for  analyzing  interlaminar  force-deformation  behavior  in  the  vicinity 
of  defects.  It  is  interesting  to  note  that  more  accurate  results  can  be 
obtained  on  reanalysis  by  simply  changing  the  element  specification  to  CCC  as 
in  the  P-Version  of  the  finite  element  method  [14]. 

Analysis  of  the  thick  laminate  impacted  off-center  by  a  small  sphere 
indicates  the  maximum  transverse  shear  occurs  below  the  surface  in  the  second 
ply  for  the  [±45/02/±45/02/±45/0/90]2s  layup.  This  maximum  occurs  at  the 
contact  radius  and  both  transverse  shears  in  material  coordinates  have  a  sinu¬ 
soidal  variation  around  the  perimeter.  The  fiber  stresses  also  vary  sinusoid¬ 
ally,  but  with  one-half  the  period  of  the  transverse  shears.  Only  the  in¬ 
plane  shear  stress  shows  a  noticeable  effect  of  the  impact  site  being  off- 
center.  The  edge  closest  the  spar  boundary  has  slightly  higher  ply  shear. 

This  also  occurs  in  the  ply  shear  strains  with  the  relative  increase  being 
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on  the  order  of  10  percent.  These  results  were  obtained  after  several  analysis 
failures  caused  by  mixing  CCL  and  CCC  elements  in  the  thickness  direction. 

Based  on  these  results  and  those  from  the  disbonded  laminate  the  following 
additional  analyses  are  recommended: 

•  Analyze  the  disbonded  laminate  without  any  symmetry 
assumptions. 

•  Analyze  the  disbonded  laminate  with  twice  as  many 
plies,  [0/+45/02/±45/0]2c,  with  the  disbond  in  the 
same  location. 

•  Analyze  the  impact  site  with  the  focus  on  the  bottom 
ply  group. 

Now  that  accurate  modeling  techniques  for  finite  element  analysis  of 
interlaminar  behavior  have  been  developed,  a  combined  test/analysis  study  of 
the  effects  of  defects  is  recommended.  It  would  be  possible  to  evaluate  energy 
release  rate  and  possibly  other  parameters  as  a  measure  of  how  large  a  delami¬ 
nation  can  become  before  unstable  propagation  takes  place.  Pre-test  analyses 
could  be  used  to  define  the  critical  load  conditions  and  to  help  locate  instru¬ 
mentation.  In  order  to  facilitate  this  work  it  would  be  desirable  to  interface 
the  PATCHES-III  program  with  a  composite  material  synthesis  program  [15],  [16] 
for  pre-  and  postprocessing  of  the  finite  element  results.  This  would  allow 
basic  ply  properties  and  layup  data  as  input  and  automate  the  recovery  of  ply 
stress-strain  data  from  composite  stress-strain  results.  It  would  also  allow 
Tsai-Wu,  Hill  and  other  failure  criteria  to  be  applied  to  the  PATCHES-III 
results. 
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APPENDIX 

The  PATCHES-III  input  data  for  the  two  laminates  analyzed  in  this 
report  and  data  modifications  necessary  to  analyze  three  additional  impact 
cases,  Table  A-l,  are  provided  in  this  appendix.  The  structural  response 
model  that  provides  displacement  boundary  conditions  for  the  impact  site  model 
is  included  for  each  case.  Also,  a  schematic  of  a  more  detailed  model  of  the 
impact  site  that  includes  the  entire  12  plies  in  the  basic  repeating  ply  group 
is  provided. 

TABLE  A-l.  PATCHES-III  LAMINATE  IMPACT  CASES 


Case 

Impact 

Site 

Sphere 

Radius 

Contact 

Radius 

Status 

1 

Z1 

=  1.91  cm 

0.635  cm 

1 .12  mm 

Analyzed 

2 

Z1 

=  1.91  cm 

2.540  cm 

1.78  mm 

Modeled 

3 

Z1 

=  10.16  cm 

0.635  cm 

1 .12  mm 

Modeled 

4 

Z1 

=  10.16  cm 

2 . 540  cm 

1 .78  mm 

Modeled 
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PATCHES- III  INPUT  DATA 
DISBONDED  LAMINATE  MODEL 


i 

I 


A- 2 


HP?PAT 
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PATCHES-III  INPUT  DATA 
LAMINATE  IMPACT 
STRUCTURAL  RESPONSE  MODEL 


Case  1.  (See  Computer  Listing) 

Case  2.  Change  gridpoints  9,  11,  13,  14,  17,  19,  21,  23. 

GRID,  9,  ,-0.356-2,-0.356-2 

GRID, 11,  ,  0.356-2,-0.356-2 

GRID, 13,  ,-0.178-2,-0.178-2 

GRID, 15,  ,  0.178-2,-0.178-2 

GRID, 17,  ,-0.178-2,  0.178-2 

GRID, 19,  ,  0.178-2,  0.178-2 

GRID, 21,  ,-0.356-2,  0.356-2 

GRID, 23,  ,  0.356-2,  0.356-2 

Case  3.  Change  gridpoints  5,  7,  29,  31  and  redefine  element  number  11  data. 


GRID,  5 

10 

.16-2,- 

-2.0-2 

GRID,  7 

-10 

.16-2,- 

-2.0-2 

GRID, 29 

10 

.16-2, 

2.0-2 

GRID, 31 

-10 

.16-2, 

2.0-2 

PATCHQ 

,11 

7, 

31,25 

1 

CPDE3 

,11 

7, 

31,25 

1 

+C11  , 

CCL 

8, 

32,26 

2 

SDC10  , 

10 

11, 

123,  7 

8,32, 

,31 

SDC10  , 

10 

11, 

2,  7 

1,  2, 

,  8 

SDC10  , 

10 

11, 

2,31 

25,26, 

,32 

SDC10  , 

10 

10, 

123,  5 

6,30, 

,29 

(delete)  SDC10  , 

10 

8, 

123,  1 

2,26, 

,25 

Case  4.  Combine  the  changes  for  Cases  2  and  3. 
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LAMINATE  STRUCTURAL  RESPONSE  MODEL 


FOR  CASES  3  AND  4 
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PATCHES- III  INPUT  DATA 
LAMINATE  IMPACT 
IMPACT  SITE  MODEL 

Case  1.  (See  Computer  Listing) 

Case  2.  (See  Computer  Listing  Page  A-29  for  gridpoint  changes). 
Case  3.  Same  as  Case  1  with  new  structural  response  input  data. 
Case  4.  Same  as  Case  2  with  new  structural  response  input  data. 
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PATCHES- III  SCHEMATIC 
LAMINATE  IMPACT 
IMPACT  SITE  12  PLY  MODEL 

The  input  data  for  this  model  are  a  direct  extension  of  the  4  ply 
model  with  32  finite  elements  increased  to  88  finite  elements.  If  the  find¬ 
ings  of  the  present  report  are  used,  then  one  plane  of  symmetry  (e^e^)  can 
be  used  to  reduce  the  number  of  elements.  In  this  instance  the  topologically 
equivalent  model  would  have  66  finite  elements.  To  reduce  this  number  further 
will  require  a  change  to  the  basic  layout  used  in  this  report. 
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(Attn:  Mr.  A.  Gustafson) .  1 

USAMATRESAG,  Watertown,  MA 

(Attn:  Dr.  E.  Lenoe) .  1 

USARESOFC,  Durham,  NC  27701  .  1 

Non-Government  Agencies 

Avco  Aero  Structures  Division,  Nashville,  TN  37202 

(Attn:  Mr.  W.  Ottenville)  .  1 

Batelle  Columbus  Laboratories,  Metals  and  Ceramics 

Information  Center,  505  King  Avenue,  OH  43201  .  1 

Bell  Aerospace  Company,  Buffalo,  NY  14240 

(Attn:  Zone  1-85,  Mr.  F.  M.  Anthony) .  1 

Bell  Helicopter  Company,  Fort  Worth,  TX  76100 

(Attn:  Mr.  Charles  Harvey) .  1 

Bendix  Products  Aerospace  Division,  South  Bend,  IN  46619 

(Attn:  Mr.  R.  V.  Cervelli) .  1 

Boeing  Aerospace  Company,  P.0.  Box  3999,  Seattle,  WA  98124 

(Attn:  Code  206,  Mr.  R.  E.  Horton) .  1 

Boeing  Company,  Renton,  Washington  98055 

(Attn:  Dr.  R.  June) .  1 

Boeing  Company,  Vertol  Division,  Philadelphia,  PA  19142 

(Attn:  Mr.  R.  L.  Pinckney,  Mr.  D.  Hoffstedt) .  2 

Boeing  Company,  Wichita,  KS  67210 

(Attn:  V.  Reneau/MS  16-39) .  1 

Cabot  Corp.,  Billerica  Research  Center,  Billerica,  MA  01821.  ...  1 

Drexel  University,  Philadelphia,  PA  19104 

(Attn:  Dr.  P.  C.  Chou) .  1 

(Attn:  Dr.  A.  S.  D.  Wang) .  1 

Effects  Technology,  Inc.,  5383  Hollister  Avenue,  P.0.  Box  30400, 

Santa  Barbara,  CA  93105  (Attn:  Robert  Globus) .  1 

E.  7.  DuPont  Company,  Wilmington,  DE  19898 

(Attn:  Dr.  J.  Pigoiacampi) .  1 

Fairchild  Industries,  Hagerstown,  MD  21740 

(Attr':  Mr.  D.  Buck) .  1 

Georgia  Institute  of  Technology,  Atlanta,  GA 

(Attn:  Prof.  W.  H.  Horton) .  1 

General  Dynamics/Convair,  San  Diego,  CA  92138 

(Attn:  Mr.  D.  R.  Dunbar,  W.  G.  Scheck) .  2 

General  Dynamics,  Fort  Worth,  TX  76101 

(Attn:  Mr.  J.  A.  Fant  (Mail  Zone-2844)) .  1 

(Attn:  Mr.  D.  Wilkins  (Composite  Structures  Eng.  Dept.)).  ...  1 
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General  Electric  Company,  Philadelphia,  PA  19101 

(Attn:  Mr.  L.  McCreight) .  1 

Great  Lakes  Carbon  Corporation,  New  York,  NY  10017 

(Attn:  Mr.  W.  R.  Benn,  Mgr.,  Market  Development) .  1 

School  of  Engineering  &  Applied  Science,  Materials  Research 
Laboratory,  Washington  University,  Campus  Box  1087, 

St.  Louis,  M0  63130  (Attn:  T.  Hahn)  .  1 

University  of  Delaware,  Mechanics  &  Aerospace  Eng.  Dept., 

Evans  Hall,  Newark,  DE  19711  (Attn:  Dr.  R.  B.  Pipes) .  1 

Grumman  Aerospace  Corporation,  Bethpage,  L.I.,  NY  11714 

(Attn:  Mr.  R.  Hadcock,  Mr.  S.  Dastin) .  2 

Hercules  Powder  Company,  Inc.,  Cumberland,  MD  21501 

(Attn:  Mr.  D.  Hug) .  1 

H.  I.  Thompson  Fiber  Glass  Company,  Gardena,  CA  90249 

(Attn:  Mr.  N.  Myers) .  1 

ITT  Research  Institute,  Chicago,  IL  60616 

(Attn:  Mr.  K.  Hofar) .  1 

J.  P.  Stevens  &  Co.,  Inc.,  New  York,  NY  10036 

(Attn:  Mr.  H.  I.  Shulock) .  1 

Kaman  Aircraft  Corporation,  Bloomfield,  CT  06002 

(Attn:  Tech.  Library)  .  1 

Lehigh  University,  Bethlehem,  PA  18015 

(Attn:  Dr.  G.  C.  Sih) .  1 

Lockheed-California  Company,  Burbank,  CA  91520 

(Attn:  Mr.  E.  K.  Walker,  Mr.  Vaughn) .  2 

Lockheed-Georgia  Company,  Marietta,  GA  30063 

(Attn:  Technical  Information  Dept.,  Dept.  72-34,  Zone  26)  .  .  1 

Vought  Corporation,  Dallas,  TX  75222 

(Attn:  Mr.  0.  E.  Dhonau/2-53442,  C.  R.  Foreman) .  2 

Martin  Company,  Baltimore,  MD  21203 

(Attn:  Mr.  J.  E.  Pawken) .  1 

Materials  Sciences  Corp.,  Blue  Bell,  PA  19422 .  1 

McDonnell  Douglas  Corporation,  St.  Louis,  MO  63166 

(Attn:  Mr.  Harold  Dill,  C.  Stenberg,  R.  Garret)  .  3 

McDonnell  Douglas  Corporation,  Long  Beach,  CA  90801 

(Attn:  H.  C.  Schjulderup,  G.  Lehman) .  2 

Minnesota  Mining  &  Manufacturing  Company,  St.  Paul,  MN  55104 

(Attn:  Mr.  W.  Davis) .  1 

Northrop  Aircraft  Corp.,  Norair  Division,  Hawthorne,  CA  90250 
(Attn:  Mr.  R.  D.  Hayes,  Mr.  D.  Stansbarger, 

Mr.  R.  C.  Isemann,  Mr.  R.  M.  Veritte) .  4 

Rockwell  International,  Columbus,  OH  43216 

(Attn:  Mr.  F.  Kaufman) .  2 

Rockwell  International,  Los  Angeles,  CA  90053 

(Attn:  Dr.  L.  Lackman) .  1 

Rockwell  International,  Tulsa,  OK  74151 

(Attn:  Mr.  E.  Sanders,  Mr.  J.  H.  Powell) .  2 
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Ownes  Corning  Fiberlass,  Granville,  OH  43023 

(Attn:  Mr.  D.  Mettes) .  1 

Rohr  Corporation,  Riverside,  CA  92503 

(Attn:  Dr.  F.  Riel  and  Mr.  R.  Elkin) .  2 

Ryan  Aeronautical  Company,  San  Diego,  CA  92112 

(Attn:  Mr.  R.  Long) .  1 

Sikorsky  Aircraft,  Stratford,  CT  06497 

(Attn:  Mr.  J.  Ray) .  1 

University  of  Oklahoma,  Norman,  OK  93069 

(Attn:  Dr.  G.  M.  Nordby) .  1 

Union  Carbide  Corporation,  Cleveland,  OH  44101 

(Attn:  Dr.  H.  F.  Volk) .  1 

University  of  Wyoming,  Laramie,  WY  82071 

(Attn:  Dr.  D.  F.  Adams) .  1 

Virginia  Tech,  Blacksburg,  VA  24061 

(Attn:  Dr.  K.  Reifsnider) .  1 

Compositek  Engineering  Corp. ,  6925-1  Aragon  Circle,  Buena  Park, 

CA  90620  (Attn:  Mr.  J.  V.  Noyes) .  1 

Lockheed-Cal ifornia  Co.,  Rye  Canyon  Research  Lab,  Burbank 

CA  91520  (Attn:  Mr.  Don  E.  Pettit) .  1 

Villanova  University,  Philadelphia,  PA  19085 

(Attn:  Dr.  P.  V.  McLaughlin) .  1 


1 


